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Abstract. The polynomials that arise as coefficients when a power series is raised to the 
power x include many important special cases, which have surprising properties that are 
not widely known. This paper explains how to recognize and use such properties, and it 
closes with a general result about approximating such polynomials asymptotically. 

A family of polynomials F (x), Fi(x), F 2 (x), . . . forms a convolution family if F n (x) has degree < n 
and if the convolution condition 

F n {x + y) = F n (x)F {y) + F n _ 1 {x)F 1 {y) + ■■■ + F 1 {x)F n _ l {y) + F (x)F n (y) 

holds for all x and y and for all n > 0. Many such families are known, and they appear frequently 
in applications. For example, we can let F n (x) = x n /n\; the condition 

(x + y) n _ ^ x k y n ~ k 
n\ ^ k\ (n- k)l 

is equivalent to the binomial theorem for integer exponents. Or we can let F n (x) be the binomial 
coefficient ; the corresponding identity 



.r +■ U \ ( •'' ^ / H 

k=0 



n J \kj \n — k 



is commonly called Vandermonde's convolution. 

How special is the convolution condition? Mathematica will readily find all sequences of 
polynomials that work for, say, < n < 4: 

F[n_,x_] :=Sum[f [n, j]x"j ,{j ,0,n}] /n! 

conv[n_] : =LogicalExpand [Series [F [n,x+y] ,{x,0,n},{y,0,n]-] 

==Series[Sum[F[k,x]F[n-k,y] ,{k,0,n>] ,{x,0,n},{y,0,n}]] 
Solve [Table [conv [n] , {n , , 4}] , 

[Flatten [Table [f [i , j] , {i ,0 ,4}, { j , ,4}] ] ] ] 

Mathematica replies that the i^'s are either identically zero or the coefficients of F n (x) = (f n o + 
fnix + f n2 x 2 H h f nn x n )/n\ satisfy 

foo = 1 j /io = /20 = /30 = /40 = , 

/22 = fu , fz2 = 3/n/21 ; f33 = fit i 

/42 = 4/n/3i + 3/li , / 43 = 6/^/21 , /44 = fll ■ 



1 



This allows us to choose /n, f 2 i, /31, and /41 freely. 

Suppose we weaken the requirements by asking only that the convolution condition hold when 
x = y. The definition of conv then becomes simply 

conv [n_] : =LogicalExpand [Series [F [n , 2x] , {x , , n}] 

==Series [Sum[F[k,x]F[n-k,x] ,{k,0,n}] ,{x,0,n}]] 

and we discover that the same solutions occur. In other words, the weaker requirements imply that 
the strong requirements are fulfilled as well. 

In fact, it is not difficult to discover a simple rule that characterizes all "convolution families." 

Let 

F(z) = l + F lZ + F 2 z 2 + F 3 z 3 + ■■■ 
be any power series with F(0) = 1. Then the polynomials 

F n (x) = [z n ]F(z) x 

form a convolution family. Conversely, every convolution family arises in this way or is identically 
zero. (Here the notation l [z n ] expr' stands for what Mathematica calls Coefficient [ expr, z,n] .) 

Proof. Let f(z) = InF(z) = f\z + f 2 z 2 /2\ + f 3 z 3 /3l + ■ ■ ■ . It is easy to verify that the coefficient 
of z n in F(z) x is indeed a polynomial in x of degree < n, because F(z) x = e x ^ z ^ = exp(x/iz + 
xf 2 z 2 /2\ + xf 3 z 3 / '3! H ) expands to the power series 

E ki+k 2 +k 3 +- AV2 2 /3 3 • • • k 1 +2k 2 +3k 3 + --- . 

A4!2! fc 2fc 2 !3! fc 3 k 3 l ... 

fcl,fc2,fc3,..->0 

when k\ + 2k 2 + 3/^3 + • • • = n the coefficient of z n is a polynomial in x with terms of degree 
ki + k 2 + &3 + • • • < n. This construction produces a convolution family because of the rule for 
forming coefficients of the product F{z) x+y = F{z) x F{z) v . 

Conversely, suppose the polynomials F n (x) form a convolution family. The condition F o (0) = 
F (0) 2 can hold only if F (x) = or F (x) = 1. In the former case it is easy to prove by induction 
that F n (x) = for all n. Otherwise, the condition F n (0) = 2F n (0) for n > implies that F n (0) = 
for n > 0. If we equate coefficients of x k on both sides of 

F n {2x) = F n {x)F {x) + F n _i(a;)Fi(x) + • • • + + F (x)F n {x) , 

we now find that the coefficient f n k of x k in n\F n (x) is forced to have certain values based on 
the coefficients of F\(x), . . . , F n _i(x), when k > 1, because 2 k f n k occurs on the left and 2f n k 
on the right. The coefficient f n \ can, however, be chosen freely. Any such choice must make 
F n (x) = [z n ] exp(x/nz + xf 2 \z 2 /2\ + x/ 3 iz 3 /3! + •••), by induction on n. 

Examples. The first example mentioned above, F n (x) = x n /n\, comes from the power series 
F(z) = e z ; the second example, F n (x) = (^), comes from F(z) = 1 + z. Several other power series 
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are also known to have simple coefficients when we raise them to the power x. If F(z) = 1/(1 — z), 
for instance, we find 

[z n ](i- z )-* = (~ x \-iy fj ' + "~ v 



x n l \ n 

It is convenient to use the notations 

x— = x(x — 1) ... (x — n + 1) = x\/{x — ri)\ 
x" = x(x + 1) . . . (x + n - 1) = r(x + n)/T(x) 

for falling factorial powers and rising factorial powers. Since = x—jn\ and ( a:+ ™ _1 ) = x n /n\, 
our last two examples have shown that the polynomials x—/n\ and x n /n\ form convolution families, 
corresponding to F(z) = 1 + z and F(z) = 1/(1 — z). Similarly, the polynomials 

x(x — s)(x — 2s) ... (x — (n — 1) s) 

F n{x) = j 

n! 

form a convolution family corresponding to (1 + sz) 1 ^ when s^O. 

The cases F(z) = 1 + z and F(z) = 1/(1 — z) are in fact simply the cases t = and t = 1 of a 
general formula for the binomial power series B t (z), which satisfies 

B t (z) = l + zB t (z) t . 

When t is any real or complex number, exponentiation of this series is known to yield 

,„,„,„ fx + tn\ x x(x + tn — l)...(x + tn — n + l) 

\z } BAzY = = — 7 ; 

1 i y J \ n J x + tn n\ 

see, for example, [Graham et al 1989, section 7.5, example 5], where a combinatorial proof is given. 
The special cases t = 2 and t = — 1, 



i_ v /T34j i 



B 2 (z) = ^ = 1 + z + 2z 2 + 5z 3 + 14z 4 + 42z 5 + • • • , 

B _^ z ) = 1 + V 1 + 4 ^ = ! + z _ z 2 + 2z 3 _ 5z 4 + Mz 5 _ . . . ^ 

in which the coefficients are the Catalan numbers, arise in numerous applications. For example, 
B 2 {z) is the generating function for binary trees, and B±(— z) is the reciprocal of B 2 (z). We can get 
identities in trigonometry by noting that B 2 {(^ sin#) 2 ) = sec 2 (#/2). Furthermore, if p and q are 
probabilities with p + q = 1, it turns out that B 2 (pq) = 1/ max(p, <j). The case t = 1/2, 



„ , , z + ^/IT^X z 2 z 3 z 5 2z 7 5z 9 Uz 11 
Bx/2{z)=^ j =l + z + Y + ¥ - ¥ + w - w + - rg , 

is another interesting series in which the Catalan numbers can be seen. The convolution polynomials 
in this case are the "central factorials" x(x + § — 1)— -/n\ [Riordan 1968, section 6.5], also called 
Steffensen polynomials [Roman and Rota 1978, example 6]. 
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The convolution formula corresponding to B t (z), 

x + y + tn\ x + y fx + tk\ x (y + i(n — fc)^ y 



n J x + y + tn ^— ^ \ A; J x + tk \ n — k J y + t(n — fc) 

is a rather startling generalization of Vandemonde's convolution; it is an identity for all x, y, t, 
and n. 

The limit of B t (z/ty as t — ► oo is another important function T(z)/z; here 



n-l o„3 



3z 3 8z 4 125z 5 



r « = E — — ^ 2 + - + ir + — + 

n>l 

is called the tree function because n n ~ 1 is the number of labeled, rooted trees. The tree function 
satisfies 

T(z) = ze T ^ , 
and we have the corresponding convolution family 

\ra-l 



[z n ] 

The related power series 



( T{Z) Y \z n ]c* T ^ X{X + ny 
\ z J n\ 



m// . 1 «™^ n 2 9z 3 32z 4 625z 5 

1 + ^> = T^T(T)=E— = l + + - + — + — + ... 

V ; n>0 

defines yet another convolution family of importance: We have 

r ral \ _ tn{ x ) 

[ J {l-T(z)) x n\ ' 

where t n (x) is called the tree polynomial of order n [Knuth and Pittel 1989]. The coefficients of 
t n (%) = t n \x + t n2 x 2 + • • • + t nn x n are integers with combinatorial significance; namely, t n k is the 
number of mappings of an n-element set into itself having exactly k cycles. 

A similar but simpler sequence arises from the coefficients of powers of e ze " : 



n\ [z n ]e xze * =J2(^jk n ~ k x k . 



k=0 

The coefficient of x k is the number of idempotent mappings of an n-element set into itself, having 
exactly k cycles [Harris and Schoenfeld 1967]. 

If the reader still isn't convinced that convolution families are worthy of detailed study, well, 
there's not much hope, although another example or two might clinch the argument. Consider the 
power series 

, , b n z n z 2z 2 5z 3 I5z 4 52z 5 

e =^^r = 1 + i! + ^r + ^r + T + ^r + --- ; 
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these coefficients b n are the so-called Bell numbers, the number of ways to partition sets of size n 
into subsets. For example, the five partitions that make 63 = 5 are 

{1,2,3}, {1}{2,3}, {1,2}{3}, {1,3}{2}, {1}{2}{3} . 

The corresponding convolution family is 

[Z Je n! 

where the Stirling number {^} is the number of partitions into exactly k subsets. 

Need more examples? If the coefficients of F(z) are arbitrary nonnegative numbers with a finite 
sum S, then F(z)/S defines a discrete probability distribution, and the convolution polynomial 
F n (x) is S x times the probability of obtaining the value n as the sum of x independent random 
variables having that distribution. 

A derived convolution. Every convolution family {F n (x)} satisfies another general convolution 
formula in addition to the one we began with: 



(x + y)^k F k (x) F n _ k (y) =xnF n (x + y). 



k=0 

For example, if F n (x) is the convolution family corresponding to powers of B t (z), this formula says 
that 

(x )^k( x + t ^\ x fy + t ( n ~ y _ xn ( x + y + tn \ x + y 

iho ^ k ' X + tk ^ n ~ k / y + t(n-k) ~ Xn \ n J x + y + tn' 
it looks messy, but it simplifies to another amazing identity in four parameters, 

'x + t(n — k)\ (y + tk\ y (x + y + tn 



E 

k=0 



n — k J \ k J y + tk \ n 



if we replace nbyn+1, kby n + 1 — k, and x by x — t + 1. This identity has an interesting history 
going back to Rothe in 1793 see [Gould and Kaucky 1966]. 

The alternative convolution formula is proved by differentiating the basic identity F(z) x = 
Sn>o Fn(%)z n with respect to z and multiplying by z: 

xzF'(z) F{zf- 1 = n F n{x) z n . 

n>0 

Now X^=o kFk{x) F n _fc(y) is the coefficient of z n in xzF'(z) F(z) x+y ~ 1 , while nF n (x + y) is the 
coefficient of z n in (x + y)zF'(z) F{z) x+ y- x . Q.E.D. 

Convolution and composition. Once upon a time I was trying to remember the form of a 
general convolution family, so I gave Mathematica the following command: 

Simplify [Series [ (1+Sum [A [k]z~k,{k,4>] )~x,{z,0, 4}] ] 
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The result was a surprise. Instead of presenting the coefficient of z n as a polynomial in x, Math- 
ematica chose another form: The coefficient of z 2 , for example, was ^A 2 x(x — 1) + A2X. In the 
notation of falling factorial powers, Mathematical response took the form 

1 + A x xz + {\A\x±- + A 2 x) z 2 + (\A\x^ + AiA 2 x- + A 3 x) z 3 

+ {±Ajx± + \A\A 2 x 3 + (Ax A 3 + \A 2 2 ) x± + A 4 x) z 4 + 0{zf . 

I wasn't prepared to work with factorial powers, so I tried another tack: 
Simplify [Series [Exp [Sum [a [k] z~k , {k , 4>] x] , {z , , 4}] ] 

This time I got ordinary polynomials in x, but — lo and behold — they were 

1 + a\xz + (^a 2 x 2 + a 2 x) z 2 + (^a 3 x 3 + a\a 2 x 2 + a^x) z 3 

+ (^afx 4 + \a 2 a 2 x 3 + (aia 3 + \a 2 )x 2 + a 4 x) z 4 + O(z) 5 . 

The result was exactly the same as before, but with a's in place of A's, and with normal powers in 
place of the factorials! 

So I learned a curious phenomenon: If we take any convolution family and replace each power 
x k by x—, we get another convolution family. (By the way, the replacement can be done in Math- 
ematica by saying 

Expand [F [n,x] ] / .Power [x,k_] ->k ! Binomial [x,k] ; 

expansion is necessary in case F n (x) has been factored.) 

The proof was not difficult to find, once I psyched out how Mathematica might have come up 
with its factorial-based formula: We have 

e^ = l + f(z)x + ^x 2 + ^fx 3 + ... , 

and furthermore 

(l + f(z)) x = l + f(z)x + l^x 2 - + l^fx 3 - + ... . 

Therefore if we start with the convolution family F n (x) corresponding to F(z) = e^ z \ and replace 
each x k by x-, we get the convolution family corresponding to 1 + f(z) = 1 + In F(z). 

A similar derivation shows that if we replace x k by the rising factorial power x k instead, we 
get the convolution family corresponding to l/(l — f(z)) = 1/(1 — In F(z)). In particular, if we 
begin with the family F n (x) = x(x + n) n ~ 1 /nl corresponding to T(z)/z = e T ( z \ and if we replace x k 
by x k to get 

this must be [z n ] (l — T(z)) x = t n (x)/n\, the tree polynomial. 
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Indeed, we can replace each x k by k\ Gk(x), where {Gk(x)} is any convolution family whatever! 
The previous examples, x— and x k , are merely the special cases k\ and k\ ( x+, f,~ 1 ) corresponding 
to two of the simplest and most basic families we have considered. In general we get 

1 + f(z) + ^ 2! G 2 (x) + 3! G 3 (x) + ■ ■ ■ , 

which is none other than G(/(z)) x = G(lnF(z)) x . 

For example, G k (x) = ( x+ ^ k ) -^k = x{x + 2k-l)^-/k\ is the family corresponding to B 2 {z). 
If we know the family F n (x) corresponding to we can replace x k by x(x + 2k — 1)^—5-, thereby 
obtaining the family that corresponds to B 2 (f(z)) = (l + y/l — 4/(z) ) /2f(z). 

Convolution matrices. I knew that such remarkable facts must have been discovered before, 
although they were new to me at the time. And indeed, it was not difficult to find them in books, 
once I knew what to look for. (Special cases of general theorems are not always easy to recognize, 
because any particular formula is a special case of infinitely many generalizations, almost all of 
which are false.) 

In the special case that each polynomial F n (x) has degree exactly n, i.e., when f 1 / 0, the 
polynomials n\ F n (x) are said to be of binomial type [Mullin and Rota 1970]. An extensive theory 
of such polynomial sequences has been developed [Rota et al 1973] [Garsia 1973] [Roman and Rota 
1978], based on the theory of linear operators, and the reader will find it quite interesting to compare 
the instructive treatment in those papers to the related but rather different directions explored in 
the present work. A comprehensive exposition of the operator approach appears in [Roman 1984]. 
Actually, Steffensen had defined a concept called poweroids, many years earlier [Steffensen 1941], 
and poweroids are almost exactly the same as sequences of binomial type; but Steffensen apparently 
did not realize that his poweroids satisfy the convolution property, which we can readily deduce 
(with hindsight) from equations (6) and (7) of his paper. 

Eri Jabotinsky introduced a nice way to understand the phenomena of convolution polynomials, 
by considering the infinite matrix of coefficients f n k [Jabotinsky 1947]. Let us recapitulate the 
notation that was introduced informally above: 

e xf{z) = F ^x = i + f 1 (x)z + F 2 {x) z 2 + ■■■ ; 
F„(x) = (f nl x + f n2 x 2 + ■■■ + fnnX n )/n\ ; 
f(z) = fiz + f 2 z 2 /2\ + f 3 z 3 /3\ + ■ ■ ■ . 

Then Jabotinsky 's matrix F = (f n k) is a lower triangular matrix containing the coefficients of 
n\ F n (x) in the nth row. The first few rows are 

h 

h fx 

h 3/i/ 2 /1 

h 4A/3 + 3/I 6/2/2 ft, 
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as we saw earlier. In general, 

f _ \ 11 ^| fki rk 2 fk 3 

J nk iifci fe 1 !2! fc 2 /c 2 !3! fc 3 k 3 \ ... 1 2 3 "' 

summed over all k\ , k 2 , k 3 , . . . > with 

h + k 2 + k 3 -\ = k , ki + 2k 2 + 3k 3 H = n . 

(The summation is over all partitions of the integer n into k parts, where kj of the parts are equal 
to j.) We will call such an array a convolution matrix. 

If each original coefficient fj is an integer, all entries of the corresponding convolution matrix 
will be integers, because the complicated quotient of factorials in the sum is an integer — it is the 
number of ways to partition a set of n elements into k subsets with exactly kj of the subsets having 
size j. Given the first column we can compute the other columns from left to right and from top 
to bottom by using the recurrence 

n — k-\-l / 

n — l 



fnk ~ Yl (j _ Ijfj f(n-j)(k-l) 



This recurrence is based on set partitions on which the element n occurs in a subset of size j: 
There are { n Z\) ways to choose the other j — 1 elements of the subset, and the factor /( n _j)(fc_i) 
corresponds to partitioning the remaining n — j elements into k — 1 parts. 
For example, if each fj = 1, the convolution matrix begins 



1 

1 1 
1 3 1 
17 6 1 
1 15 25 10 



These are the numbers {^} that Mathematica calls StirlingS2 [n,k] ; they arose in our example 
of Bell numbers when f(z) = e z — 1. Similarly, if each fj = (j — 1)!, the first five rows are 



1 

1 1 

2 3 1 

6 11 6 1 

24 50 35 10 1 



Mathematica calls these numbers (-1) " (n-k)StirlingSl [n,k] . In this case f(z) = ln(l/(l — z)) , 
and F n (z) = The signed numbers StirlingSl [nk] , 



1 

-1 1 

2 -3 1 

-6 11 -6 1 

24 -50 35 -10 
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correspond to f(z) = ln(l + z) and F n (z) = (*). In general if we replace z by az and x by (3x, 
the effect is to multiply row n of the matrix by a n and to multiply column k by f5 k . Thus when 
/3 = a -1 , the net effect is to multiply /„& by a n_fc . Transforming the signs by a factor (— l) n_fc 
corresponds to changing F(z) to l/F(—z) and /(z) to —f(—z). Therefore the matrix that begins 

1 

-1 1 

1 -3 1 
-1 7-6 1 

1 -15 25 -10 1 

corresponds to f(z) = 1 — e~ z . 

Let's look briefly at some of our other examples in matrix form. When F(z) = Bt(z), we have 
fj = (tj — 1)^-, which is an integer when t is an integer. In particular, the Catalan case t = 2 
produces a matrix that begins 

1 

3 1 

20 9 1 

210 107 18 1 

3024 1650 335 30 1 . 
When t = 1/2, we can remain in an all- integer realm by replacing z by 2z and x by x/2. Then 

fj = when j is even, while /iy+i = (— 1) J (2j — l)!! 2 : 

1 

1 
-1 1 
0-4 1 ; 
9 -10 1 

If we now replace z by iz and x by xji to eliminate the minus signs, we find that f(z) = arcsin z, 
because ln(iz + y/l — z 2 ) = id when z = sm.6. Thus we can deduce a closed form for the coefficients 
of e xaTCsinz = B 1/2 {2iz) x /^: 

n\ [z n ] e ™ inz = ( 2i )"-i x (?L + " _ A . . . (* _ ^ + A 

Ji 2 (x 2 + 2 2 ) ... (x 2 + (n-2) 2 ), n even; 

" { x(x 2 + l 2 )(x 2 + 3 2 ) . . . (x 2 + ((n - 2) 2 )) , n odd. 
This remarkable formula is equivalent to the theorem of [Gomes Teixeira 1896]. 

If fj = 2 1-7 when j is odd but fj = when j is even, we get the convolution matrix corre- 
sponding to e 2x smh ( z / 2 ) : 



1 

1 



1 



1 

4 

10 1 

— - 1 

16 2 

1 5 1 

J_ n 91 — 1 

64 u 16 U 4 u 1 

1 21 14 
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Again we could stay in an all-integer realm if we replaced z by 2z and x by x/2; but the surprising 
thing in this case is that the entries in even-numbered rows and columns are all integers before we 
make any such replacement. The reason is that the entries satisfy f n k = fc 2 /(n-2)fc/4 + /( n -2)(fc-2)- 
(See [Riordan 1968, pages 213-217], where the notation T(n,k) is used for these "central factorial 
numbers" f nk .) 

We can complete our listing of noteworthy examples by setting fj = Y^k=i n n ~ k ~ 1 n-; then 
we get the coefficients of the tree polynomials: 



1 








3 


1 






17 


9 


1 




142 


95 


18 


1 


1569 


1220 


305 


30 



The sum of the entries in row n is n™. 

Composition and iteration. Jabotinski's main reason for defining things as he did was his 

observation that the product of convolution matrices is a convolution matrix. Indeed, if F and G are 
the convolution matrices corresponding to the functions e x ^ z " > and e X9<yZ ^ we have the vector/matrix 
identities 

e x/(z) _ j = ( Z:Z 2/ 2 \,z 3 /3\,...)F(x,x 2 ,x 3 ,...) T 
e xg{z) _ 1 = ( Z:Z 2/ 2 \,z 3 /3\,...)G(x,x 2 ,x 3 ,...) T 
If we now replace x k in e x ^ z ^ by k\ Gk(x), as in our earlier discussion, we get 

(z, z 2 /2\, z 3 /3\, ...)F {G 1 (x),2\ G(x),3\ G 3 (x),...) T 

= {z, z 2 /2\, z 3 /3l, ...)FG (x,x 2 ,x 3 ,... ) T 

=(f(z), f(z) 2 /2\, f(z) 3 /3\, ...)G(x,x 2 ,x 3 ,...) T 
=e xg(f(z)) _ j _ 

Multiplication of convolution matrices corresponds to composition of the functions in the exponent. 

Why did the function corresponding to FG turn out to be g(f(z)) instead of f(g(z)) ? Jabotin- 
sky, in fact, defined his matrices as the transposes of those given here. The rows of his (up- 
per triangular) matrices were the power series f(z) k , while the columns were the polynomials 
F n (x) = [z n ] e x f( z ^; with those conventions the product of his matrices F T G T corresponded to 
f{g(z)}. (In fact, he defined a considerably more general representation, in which the matrix F 
could be U~ l FU for any nonsingular matrix U.) However, when our interest is focussed on the 
polynomials n\ F n (x), as when we study Stirling numbers or tree polynomials or the Stirling poly- 
nomials to be discussed below, it is more natural to work with lower triangular matrices and to 
insert factorial coefficients, as Comtet did [Comtet 1970, section 3.7]. The two conventions are iso- 
morphic. Without the factorials, convolution matrices are sometimes called renewal arrays [Rogers 
1978]. We would get a non-reversed order if we had been accustomed to using postfix notation 
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(z)f for functions, as we do for operations such as squaring or taking transposes or factorials; then 
g{f(z)) would be ((z)f)g. 

Recall that the Stirling numbers {^} correspond to f(z) = e z — 1, and the Stirling numbers 
[£] correspond to g(z) = ln(l/(l — z)). Therefore if we multiply Stirling's triangles we get the 
convolution matrix 

1 

2 1 
FG =661 

26 36 12 1 
150 250 120 20 1 

which corresponds to g(f(z)) = ln(l/(2 — e z )). Voila! These convolution polynomials represent 
the coefficients of (2 — e z )~ x . [Cayley 1859] showed that (2 — e z )~ 1 is the exponential generating 
function for the sequence 1, 3, 13, 75, 541, . . . , which counts preferential arrangements of n objects, 
i.e., different outcomes of sorting when equality is possible as well as inequality. The coefficient 
(fg)nk is the number of preferential arrangements in which the "current minimum" changes k times 
when we examine the elements one by one in some fixed order. (See [Graham et al 1989, exercise 
7.44].) 

Similarly, the reverse matrix product yields the so-called Lah numbers [Lah 1955], 

1 

2 1 

GF =661 ; 
24 36 12 1 
120 240 120 20 1 

here fj = j\ and the rows represent the coefficients of exp(x/(p(z))) = exp(xz+xz 2 +xz 3 + - ■ • ). In- 
deed, the convolution polynomials in this case are the generalized Laguerre polynomials Ln ^(— x), 
which Mathematica calls LaguerreL [n, -1 , -x] . These polynomials can also be expressed as 
L n (—x) — L n _i(— x); or as LaguerreL [n,-x] -LaguerreL [n-l,-x] if we say 

Unprotect [LaguerreL] ; LaguerreL [- 1 , x_] :=0; Protect [LaguerreL] 

first. The row sums 1, 3, 13, 73, 501, .. . of GF enumerate "sets of lists" [Motzkin 1971]; the coeffi- 
cients are {GF) nk = n\ [z n ] f{g{z)) k /k\ = (™) (£lj)(n - k)\ [Riordan 1968, exercise 5.7]. 

Since convolution matrices are closed under multiplication, they are also closed under exponen- 
tiation, i.e., under taking of powers. The c/th power F q of a convolution matrix then corresponds 
to g-fold iteration of the function \nF = f. Let us denote f{f(z)) by f^ 2 \z); in general, the qth. 
iterate (z) is defined to be (z)) , where (z) = z. This is Mathematical Nest[f,z,q]. 

The qth iterate can be obtained by doing 0(log q) matrix multiplications, but in the interesting 
case f'(0) = fi = 1 we can also compute the coefficients of f^(z) by using formulas in which q 
is simply a numerical parameter. Namely, as suggested by [Jabotinsky 1947], we can express the 
matrix power F q as 

(l+(F-I)) q = I+{^(F-I) + (^(F-ir + (^(F-lf + ... . 
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This infinite series converges, because the entry in row n and column k of (F — I) J is zero for all 
j > n — k. When q is any positive integer, the result defined in this way is a convolution matrix. 
Furthermore, the matrix entries are all polynomials in q. Therefore the matrix obtained by this 
infinite series is a convolution matrix for all values of q. 

Another formula for the entries of F q was presented in [Jabotinsky 1963]. Let be the 
element in row n and column k; then 
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for any m>n-k. Indeed, we have p(q) = Y^JLo P(J) (f) ( 9 m -/) ( 
of degree < m; this is a special case of Lagrange interpolation. 



- l) m J whenever p is a polynomial 



It is interesting to set q 
matrices. We have 



1/2 and compute convolution square roots of the Stirling number 
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The function z+z 2 /4+z 3 /48+z 5 /3840-7z 6 /92160+- • • therefore lies "halfway" beween z and e z -l = 
z + z 2 /2\ + z 3 /3\ + • • • , and the function z + z 2 /A + 5z 3 /48 + 5z 4 /96 + 109z 5 /3840 + 497z 6 /30720 + ■ ■ ■ 
lies halfway between z and In 1/(1 — z) = z + z 2 /2 + z 3 /3 + • • • . These half-iterates are unfamiliar 
functions; but it is not difficult to prove that z/(l - z/2) = z + z 2 /2 + z 3 /A + ■■■ is halfway 
between z and z/(l — z) = z + z 2 + z 3 + • • • . In general when f(z) = z/(l — cz k ) x l k we have 
/[«](*) = z/{l-qcz k ) l / k . 

It seems natural to conjecture that the coefficients of f^(z) are positive for q > when 
f(z) = In 1/(1 — z); but this conjecture turns out to be false, because Mathematica reports that 
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[z s ] f[<A( z ) = -llg/241920 + 0{q 2 ). Is there a simple necessary and sufficient condition on / that 
characterizes when all coefficients of are nonnegative for nonnegative q? This will happen if 
and only if the entries in the first column of 

In F = (F — I) — \{F - I) 2 + \{F -If 

are nonnegative. (See [Kuczma 1968] for iteration theory and an extensive bibliography.) 

Reversion. The case q = — 1 of iteration is often called reversion of scries, although Mathematica 
uses the more proper name InverseSeries. Given f(z) = f\z + f2Z 2 /2\ + • • • , we seek g(z) = 
f^~ l \z) such that (/(/(z)) = z. This is clearly equivalent to finding the first column of the inverse 
of the convolution matrix. 

The inverse does not exist when /i = 0, because the diagonal of F is zero in that case. 
Otherwise we can assume that /i = 1, because fig(f{z/fi)) = z when g reverts the power series 
f(z/h). 

When fi = 1 we can obtain the inverse by setting q = — 1 in our general formula for iteration. 
But Lagrange's celebrated inversion theorem for power series tells us that there is another, more 
informative, way to compute the function g = Let us set F(z) = f(z)/z = 1 + f2 z /2\ + 

f3Z 2 /3\ + • • • . Then Lagrange's theorem states that the elements of the matrix G = F~ l are 

(n- 1)! ~ , , 

where F n (x) denotes the convolution family corresponding to F{z). 

There is a surprisingly simple way to prove Lagrange's theorem, using our knowledge of con- 
volution families. Note first that 

/„* = n\ [z n x k ] = y [z n ] f(z) k = y [z n ~ k ] F(z) k ; 

therefore 

fnk = F n -k(k) ■ 

Now we need only verify that the matrix product GF is the identity, by computing its element in 
row n and column m: 

9nkfkm = Yl (k-l)l Fn ~ k (~ n ) F k-m(™) ■ 

k=m k=m 

When m = n the sum is obviously 1. When m = n — p for p > it is (n — l)!/(n — p)\ times 

n V 

k Fn-k(-n) F k _ n+p (n - p) = ^ (n - k) F k (-n) F p _ k (n - p) 

k=n—p k=0 

V V 



= n Y Fk{-n) F P - k (n - p) - ^ k F k {-n) F p - k {n - p) 

k=0 k=0 

= nF p (-p)-nF p (-p)=0 
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by the original convolution formula and the one we derived from it. The proof is complete. 



Extending the matrix. The simple formula for f nk that we used to prove Lagrange's theorem 
when fi = 1 can be written in another suggestive form, if we replace k by n — k: 

fn(n-k) = n-F k {n - k) . 

For every fixed k, this is a polynomial in n, of degree < 2k. Therefore we can define the quantity 
for all real or complex y to be y-F k (y — k); and in particular we can define f nk in this 
manner for all integers n and k, letting f nk = when k > n. For example, in the case of Stirling 
numbers this analysis establishes the well-known fact that { y l k } and [ y ^_ k ] are polynomials in y 
of degree 2k, and that these polynomials are multiples of y^- = y(y — 1) ... (y — k) when k > 0. 

The two flavors of Stirling numbers are related in two important ways. First, their matrices 
are inverse to each other if we attach the signs (— l) n_fc to the elements in one matrix: 
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This follows since the numbers {^} correspond to f(z) = e z — 1 and the numbers [^] (— 1)™ k 
correspond to g(z) = ln(l + z), as mentioned earlier, and we have g(f(z)^ = z. 
The other important relationship beween {^} and [^] is the striking identity 
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ft}- 


—n 



which holds for all integers n and k when we use the polynomial extension method. We can prove 
in fact, that the analogous relation 



fnk = (-l) fc n 9(-k)(-n) 

holds in the extended matrices F and G that correspond to any pair of inverse functions g(f{z)) = z, 
when /'(0) = 1. For we have 



n! 



(-I)"" V fc)( _ n) = (-l) n - k (-k -l)(-k-2) ... (-n) F n . k (k) = - F n . k (k) = f nk 

in the formulas above. (The interesting history of the identity {^} = [Z*] is traced in [Knuth 
1992]. The fact that the analogous formula holds in any convolution matrix was pointed out by Ira 
Gessel after he had read a draft of that paper. See also [Jabotinski 1953]; [Carlitz 1978]; [Roman 
and Rota 1978, section 10].) 

Suppose we denote the Lah numbers (fc)(fcZi)( n — k)\ by |£|. The extended matrix in that 
case has a pleasantly symmetrical property 
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—n 
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because the corresponding function f(z) = z/(l — z) satisfies /(—/(— 2)) = z. (Compare [Mullin 
and Rota 1969, section 9].) Near the origin n = k = 0, the nonzero entries look like this: 

... 1 
... 12 1 
... 36 6 1 
... 24 6 2 1 

1 

1 

2 1 
6 6 1 
24 36 12 1 



Still more convolutions. Our proof of Lagrange's theorem yields yet another corollary. Suppose 
g(f(z)) = z and f'(0) = 1, and let F(z) = f(z)/z, G(z) = g(z)/z. Then the equation 

9nk = G n - k {k) = ^ _ F n _ k (-n) 

tell us, after replacing n by n + k, that the identity 

n + k 



k 



G n (k) = F n (-n - k) 



holds for all positive integers k. Thus the polynomials G n (x) and F n (x) must be related by the 
formula 

(x + n) G n (x) = xF n {—x — n) . 

Now F n (x) is an arbitrary convolution family, and F n (— x) is another. We can conclude that if 
{F n (x)} is any convolution family, then so is the set of polynomials {xF n {x + n)/{x + n)}. Indeed, 
if F n {x) corresponds to the coefficients of F(z) x , our argument proves that the coefficients of G(z) x 
are x F n (x + n)/{x + n), where zG(z) is the inverse of the power series z/F(z): 

G(z) = F(zG(z)) , G(z/F(z))=F(z). 

The case F(z) = 1 + z and G(z) = 1/(1 — z) provides a simple example, where we know that 
F n (x) = □ and G n (x) = i^ 1 ) = xF n (x + n)/(x + n). 

A more interesting example arises when F{z) = ze z /(e z — 1) = z + z/(e z — 1) = 1 + z/2 + 
B2Z 2 /2\+B^z 4 /4l + - • • ; then F(—z) is the exponential generating function for the Bernoulli numbers. 
The convolution family for F(z) x is F n (x) = xo n (x), where cr n (x) is called a Stirling polynomial. 
(Actually cjo(x) = 1/x, but a n (x) is a genuine polynomial when n > 1.) The function G such 
that G{z/F(z)) = F(z) is G(z) = z _1 ln(l/(l — z)); therefore the convolution family for G(z) x is 
G n (x) = xF n (x + n)/(x + n) = i(j n (x + n) . 
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In this example the convolution family for e xzG ( z ) = (\ — z ) x is 



therefore 
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(n — /c) <Tfc (n) = n(n — 1) ... {n — k)a k (n) . 



We also have 



{n - /c} 


k — n 




—n 



(n — k)\ (n — k)\ 

(k — n) (k — 1 — n) ... (— n) <7fc(fc — n) . 

These formulas, which are polynomials in n of degree 2/c for every fixed k, explain why the a func- 
tions are called Stirling polynomials. Notice that a n {\) = (—l) n B n /nl; it can also be shown that 
a n (0) = -B n /(n-n\). 

The process of going from F n (x) to xF n (x + n)/(x + n) can be iterated: Another replacement 
gives xF n (x + 2n)/(x + 2n), and after t iterations we discover that the polynomials xF n (x + tn) / '(x + 
tn) also form a convolution family. This holds for all positive integers t, and the convolution 
condition is expressible as a set of polynomial relations in t; therefore xF n (x + tn)/(x + tn) is a 
convolution family for all complex numbers t. If F n (x) = [z n ] F{z) x , then xF n (x + tn)/(x + tn) = 
[z n ] T t {z) x , where Tt{z) is defined implicitly by the equation 

T t {z) =F{zT t {zf). 

In particular, we could have deduced the convolution properties of the coefficients of B t {z) x in this 
way. 

Let us restate what we have just proved, combining it with the "derived convolution formula" 
obtained earlier: 

Theorem. Let F n (x) be any family of polynomials in x such that F n (x) has degree < n. If 

n 

F n (2x) = Y,Fk(x)F n _ k (x) 

k=0 

holds for all n and x, then the following identities hold for all n, x, y, and t: 

(x + y) F n (x + y + tn) xF k{ x + tk ) y F n - k (y + t(n - k)) _ 



x + y + tn 



fc=0 



x + tk 



y + t(n — k) 



n 



F n {x + y + tn) = ^ kF k (x + tk) y F n - k (y + t(n - k)) 



x + y + tn 



k=i 



x + tk 



y + t(n — k) 



Additional constructions. We have considered several ways to create new convolution families 
from given ones, by multiplication or exponentiation of the associated convolution matrices, or by 



16 



replacing F n (x) by x F n (x + tn)/{x + tn). It is also clear that the polynomials a n F n (f3x) form a 
convolution family whenever the polynomials F n (x) do. 

One further operation deserves to be mentioned: If F n (x) and G n (x) are convolution families, 
then so is the family H n (x) defined by 

n 

H n (x) = ^2 F k (x) G n _ k {x) . 

k=0 

This is obvious, since H n (x) = [z n ] F{z) x G{z) x . The corresponding operation on matrices F = 
ifnk), G = (g nk ), H = (h nk ) is 

Kk = ^2 ( J fjig(n-j)(k-i) ■ 

If we denote this binary operation by H = F o G, it is interesting to observe that the associative 
law holds: {E o F) o G = E o (F o G) is true for all matrices E, F, G, not just for convolution 
matrices. A convolution matrix is characterized by the special property F o F = F diag(2, 4,8,...). 
The construction just mentioned is merely a special case of the one-parameter family 

n 

H%\x) = J2 F k(x) G n . k {x + tk). 

k=0 

Again, {Hn\x)} turns out to be a convolution family, for arbitrary t: We have 

n 

Y,H^{x)z n = J2 F k( x ) G n-k(x + tk)z n = Y, F k(x)G n (x + tk)z n+k 

k=0 n>k>0 n,k>0 

= Y^F k {x)z k G{z) x+tk = G{z) x F{zG{z) t ) x 

k>0 

so H n (x) = [z n ] {G(z)F{zG{z) t )) x . 

Applications. What's the use of all this? Well, we have shown that many interesting convolution 
families exist, and that we can deduce nonobvious facts with comparatively little effort once we 
know that we're dealing with a convolution family 

One moral to be drawn is therefore the following. Whenever you encounter a triangular pattern 
of numbers that you haven't seen before, check to see if the first three rows have the form 

a 

b a 2 

c Sab a 3 

for some a,b,c. (You may have to multiply or divide the nth row by n\ first, and/or reflect its 
entries left to right.) If so, and if the problem you are investigating is mathematically "clean," 
chances are good that the fourth row will look like 

d 4ac + 3b 2 6a 2 b a 4 . 
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And if so, chances are excellent that you are dealing with a convolution family. And if so, you may 
well be able to solve your problem. 

In fact, exactly that scenario has helped the author on several occasions. 



Asymptotics. Once you have identified a convolution family F n (x), you may well want to know 
the approximate value of F n (x) when n and x are large. The remainder of this paper discusses 
a remarkable general power series expansion, discovered with the help of Mathematica, which 
accounts for the behavior of F n (x) when n/x stays bounded and reasonably small as x — > oo, 
although n may also vary as a function of x. We will assume that F n (x) is the coefficient of z n in 
F(z) x , where F(0) = F'(0) = 1. 

Our starting point is the classical "saddle point method," which shows that in many cases the 
coefficient of z n in a power series P(z) can be approximated by considering the value of P at a point 
where the derivative of P(z)/z n is zero. (See [Good 1957].) In our case we have P(z) = e x ^ z \ 
where f(z) = In F(z) = z + /2-z 2 /2! + • • • ; and the derivative is zero when x f'(z) = n/z. Let this 
saddle point occur at z = s; thus, we have 

s f'(s) = n/x . 

Near s we have f(z) = f(s) + (z — s)f'(s) + 0({z — s) 2 ); so we will base our approximation on the 

assumption that the 0[(z — s) 2 ) contribution is zero. The approximation to F n (x) will be F n (x), 

where _ 

F n (x) = [z n ] exp (x f(s) +x(z-s) /'(*)) 

p x(f(s)-sf' (s)) F(s) x / n \ n 

= ~ j * n f'{sY = (-) • 

n! n! \es/ 

First let's look at some examples; later we will show that the ratio F n (x)/F n (x) is well behaved 
as a formal power series. Throughout this discussion we will let 

y = n/x ; 

our goal, remember, is to find approximations that are valid when y is not too large, as x and 
possibly n go to oo. 

The simplest example is, of course, F(z) = e z and f(z) = z; but we needn't sneeze at it 
because it will give us some useful calibration. In this case f"(z) = 0, so our approximation will 
be exact. We have s = y, hence 



n\ e /xv l x 



F n {x) = —\-\ =-{-) = — = F n (x). 
n\ \ey J n\ V e / n! 

Next let's consider the case F(z) = T(z)/z, f(z) = T(z), when we know that F n (x) 
zix + ny- 1 /n\. In this case zT'(z) = T(z)/(l - T(z)), so we have T(s)/(l - T(s)) = y or 

T(s) = ^—, s =^—e-y/^ 

i + y 1 + 2/ 
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because T(z) = ze T ^ . Therefore 

e xy/(i+ y ) / n (l + y) \ n (x + n) 



F n (x) = 



the ratio F n (x)/F n (x) = x/(x + n) = 1/(1 + y) is indeed near 1 when y is small. 
If F(z) = 1 + z we find, similarly, s = y/(l — y) and 



n\F n {x) 



1 \ x (n(l - y) 



1 — yJ \ ey J (x — n) a 
by Stirling's approximation we also have 



n \ FJx) = . X ' ., = ~, ~ ~r (1 - y)" 1/2 (l + 0{x~ x )\ . 

Again the ratio F n (x) / F n (x) is near 1. In general if F(z) = B t (z) the saddle point s turns out to 
bey(l + (t-l)y)*" 1 /(l + ty)*, and 

(x + tn) x+tn e- n 



n\F n (x) 



(x + (t- l)n) 



x+(t-l)n 



a similar analysis shows that this approximation is quite good, for any fixed t. 
We know that 



n! \ x x z J 



and that / n ( n _fc) is always a polynomial in n of degree < 2fc. Therefore if n 2 /x — > as x — ► oo, we 
can simply use the approximation F n (x) = (x n /n!)(l + 0(n 2 /x)). But there are many applications 
where we need a good estimate of F n (x) when n 2 /x — > oo while n/x — > 0; for example, x might be 
nlogn. In such cases F„(x) is close to F n (x) but x n /«! is not. 

We can express s/y as a power series in y by inverting the power series expression sf'(s) = y: 

s/y=l- hv + (4/| - / 3 )y 2 /2 + (15/ 2 / 3 - 30/| - U)y 3 /6 + ■■■ . 

From this we can get a formal series for F n (x), 

~ x™ exp(n(g/y)(l + / 2 s/2! + / 3 g 2 /3! + •••)- n) 

n[X) ~ n\ (s/y)" 

x n (,nh , 3n 2 / 2 2 -12n/|+4n/ 3 2 3 3 



We can also use the formula 



fc + fc 2 + fc3 + " 



f n fki f k 3 

Jn(n-k) ^ 2! k * k 2 l3l k3 k 3 \ . . . 2 3 
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where the sum is over all nonnegative k 2 , £3, . . . with k 2 + 2A; 3 + • • • = k, to write 



F„(x) = ^ (1 + 



nf 2 ~ f 2 + 0(x- 1 ) 



y 



+ 3n 2 fj - Wag + Anh + 33/ 2 2 - 12/ 3 + 0(Q y2 + 0(n3y3) 



These series are not useful asymptotically unless ny = n 2 /x is small. But the approximation F n {x) 
itself is excellent, because amazing cancellations occur when we compute the ratio: 



F„(x) 2 



h , 11/1 — 4/3 2 3 _! 

— y +0(y' 3 )+0(x ). 



Theorem. When = exp(z + f 2 z 2 /2\ + /3Z 3 /3! + • • • ) and the functions F n (x) and F n (x) are 
defined as above, the ratio F n (x)/F n (x) can be written as a formal power series ^ j>o c ijU 1 ' x ~ J ' > 
where y = n/x and the coefficients Cij are polynomials in f 2 , f 3 , . . . ■ 

The derivation just given shows that we can write F n (x)/F n {x) as a formal power series of the 
form £\ j>o o,ijn l x~ 3 ' , where a^- = when i > 2j; the surprising thing is that we also have = 
whenever i > j. Therefore we can let c^- = a^ i+ jy 

To prove the theorem, we let R(z) = \ + R\z + R 2 z 2 H stand for the terms neglected in our 

approximation: 

w - ( 1 + 7ff + Slr + S¥ + -) «w • 



The coefficient of z n is 



n — 1 

n 



F n (x) = F n (x) I 1 + R lS + R 2 s 2 + ^ ^ ^ i? 3 s 3 + •••); 



(n-l)(n-2) 



so the ratio F n (x)/F n (x) is equal to 



j,fc>0 



fc>0 



k 

k - j 



where Pj = -Rfc-s fc is a certain power series in s and x. The coefficients i?^ are themselves 

power series in s and x, because we have 



+ x(z — s) 



3 rw 



3! 



We know from the discussion above that 

k _. =k(k-l) ... (k-tia^k) 
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is a polynomial in k. Therefore we can write 

P = \ * 

3 U - 3. 



R(z) 



where # is the operator that takes z k i— > k z k for all k; i.e., ■dG(z) = zG'(z) for all power series 
G(z). The theorem will be proved if we can show that Pj/n^ is a formal power series in y and 
and if the sum of these formal power series over all j is also such a series. 

Consider, for example, the simplest case Po = R(s); obviously Po = 1- The next simplest case 
is P 1 = J R(z) = - l)-R(z) | z=g . It is easy to see that 

•dL = z j D j , 

where D is the differentiation operator D G(z) = G'(z), because z J D j takes z k into k 3 -z k . Therefore 

P 1 = ±s 2 R"(s) = \xs 2 f"{s). 
It follows that P\jn = ^(s/y)sf"(s) is a power series in y; it begins |/22/ + \{fz — f^V 2 + " " " ■ 



Now let's consider Pj in general. We will use the fact that the Stirling numbers 



represented in the form 



-k—ji 



can be 



k 

k - j 



Pji 



k 



+ Pj2 



k 

3 + 2 



+ ---+Pjj 



vi + h 

where the coefficients pji are the positive integers in the following triangular array: 

1 

2 3 

6 20 15 
24 130 210 105 
120 924 2380 2520 945 . 

(This array is clearly not a convolution matrix; but the theory developed above implies that the 
numbers jlpji/(i + j)! , namely 



1/2 
2/3 
3/2 
24/5 
20 



1/4 
1 

13/3 
22 



1/8 
1 1/16 
85/12 5/6 1/32 



do form the convolution matrix for the powers of exp(z/2 + z 2 /3 + z 3 /4 + ■ ■ ■ ). The expression 
[k- ] = Yli=iPji( +i) was independently discovered by [Appell 1880], [Jordan 1933], and [Ward 
1934]. The number of permutations of i + j elements having no fixed points and exactly i cycles 
is pji, an "associated Stirling number of the first kind" [Riordan 1958, section 4.4] [Comtet 1970, 
exercise 6.7].) It follows that 



. , R^+^fs) 

P-=Vl S J + 1 — — + V OS 

r 3 Pjl b (j + l)' Pd 



J+ 2 



(3 + 2)1 



+ ---+PjjS 



2j 
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Now R(z) is a sum of terms having the form 

aux\z - s) 1 , 

where I > 2i and where an is a power series in s. Such a term contributes aux % s p^i_^ to Py, so 
it contributes a^s/y) s l ~i x l ~ipj(i_j) to Pj/n j . This contribution is nonzero only if j < I < 2j. 
Since I > 2i, we have i < j; so Pjjn? is a power series in y and x~ l . 

For a fixed value of j — i, the smallest power of y that can occur in Pjjn 1 is y 2l ~i = y^~ 2 ^~ l \ 
Therefore only a finite number of terms of Y^j Pj/(— n) J contribute to any given power of y and 
This completes the proof. 

A careful analysis of the proof, and a bit of Mathematica hacking, yields the more precise 
result 

where ^4 = ^s 3 ^ -1 ^ _ f — \s 2 d2<lz — ^s 3 dl + |yd3 + ^s 3 d 2 d4 i + |sy<i4 and dk = f^(s). 
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